	! Molecular dynamics routines
	! Unit of distance = Angstroms
	! Unit of time = ps
	! Working unit of energy is Joules
	! Assumed input unit of force = eV / Angstrom
	! Output unit

	module mdroutines
	implicit none
	! Constants
	real*8, parameter :: pi = 3.14159265358979d0
	real*8, parameter :: kb_jmolk = 8.314472d0	! Ideal gas constant, J mol-1 K-1
	real*8, parameter :: kb_10jmolk = 0.8314472d0	! Ideal gas constant, 10 J mol-1 K-1
	real*8, parameter :: kb_evk = 8.617343d-5	! Boltzmann constant, eV K-1  (8.314472 / 1000.0 / 96.48536)
	! Number of particles in system
	integer :: md_natoms = 0
	! Particle masses, positions, velocities, and forces
	real*8, allocatable :: md_mass(:), md_r(:,:), md_rlast(:,:), md_v(:,:), md_vlast(:,:), md_f(:,:), md_temp(:,:)
	real*8 :: md_sumv(3)
	! Particle names
	character*8, allocatable :: md_atomnames(:)
	! Unit cell
	real*8 :: md_cell(9), md_icell(9)
	! Control
	integer :: md_nsteps = 1, md_step = 0
	real*8 :: md_timestep = 1.0d-5, md_Ttarget = 300.0d0, md_Tcalc, md_pe, md_ke, md_totaltime = 0.0d0
	contains

	subroutine md_createempty(natoms)
	implicit none
	integer :: natoms
	md_natoms = natoms
	if (md_natoms.lt.1) stop "Must call md_createarrays with a useful number of atoms."
	allocate(md_mass(md_natoms), md_atomnames(md_natoms))
	allocate(md_r(md_natoms, 3), md_v(md_natoms, 3), md_f(md_natoms, 3))
	allocate(md_rlast(md_natoms, 3), md_vlast(md_natoms, 3), md_temp(md_natoms, 3))
	! Set some initial values
	md_mass = 1.0d0
	end subroutine 

	subroutine md_createfcc(alpha,kx,ky,kz,massi,massj,massk,massl,quiet)
	use cellutils
	implicit none
	real*8 :: alpha, massi, massj, massk, massl
	integer :: kx, ky, kz, n, i, j, k
	logical :: quiet
	if (.not.quiet) write(6,"(/,a50,f8.4,a,3i3)") "Generating initial FCC lattice with alpha = ", alpha," and k(xyz) = ",kx,ky,kz
	! Do we need to reallocate arrays?
	if (md_natoms.ne.(4*kx*ky*kz)) call md_createempty(4*kx*ky*kz)
	n = 1
	do i = 0,kx-1
	  do j = 0,ky-1
	    do k = 0,kz-1
	      md_r(n,:) = (/ i*1.0d0, j*1.0d0, k*1.0d0 /)
	      md_r(n+1,:) = (/ i*1.0d0+0.5d0, j*1.0d0+0.5d0, k*1.0d0 /)
	      md_r(n+2,:) = (/ i*1.0d0+0.5d0, j*1.0d0, k*1.0d0+0.5d0 /)
	      md_r(n+3,:) = (/ i*1.0d0, j*1.0d0+0.5d0, k*1.0d0+0.5d0 /)
	      md_mass(n:n+3) = (/ massi, massj, massk, massl /)
	      n = n + 4
	    end do
	  end do
	end do
	md_r = md_r * alpha
	! Setup unit cell
	md_cell = (/ kx*1.0d0,0.0d0,0.0d0,0.0d0,ky*1.0d0,0.0d0,0.0d0,0.0d0,kz*1.0d0 /)
	md_cell = md_cell * alpha
	md_icell = md_cell
	call calc_icell(md_icell)
	! Write out some info on the generated system
	if (.not.quiet) then
	  write(6,"(a50, i5)") "Number of atoms in FCC lattice = ", md_natoms
	  write(6,"(a50,3f12.6,2(/,50x,3f12.6))") "Unit cell = ", md_cell
	  write(6,"(a50,3f12.6)") "Volume (Angstroms**3) = ", volume(md_cell)
	  write(6,"(a50,3f12.6)") "Density (g/cc) = ", (sum(md_mass)/6.02214179d23) / (volume(md_cell)/1.0d24)
	end if
	end subroutine md_createfcc

	subroutine md_freearrays()
	implicit none
	deallocate(md_mass, md_r, md_v, md_f, md_atomnames)
	deallocate(md_rlast, md_vlast, md_temp)
	end subroutine 

	subroutine md_initialise()
	use random
	implicit none
	real*8 :: t_instant, ke
	integer :: i
	! Initialise random atomic velocities according to requested temperature
	write(6,"(a,f9.2)") "Initialising particle velocities for T = ", md_Ttarget
	md_sumv = 0.0d0
	do i=1,md_natoms
	  md_v(i,1) = dexp( -(2.0d0*ran()-1.0d0) / 2.0d0 ) / dsqrt(2.0d0 * pi)
	  md_v(i,2) = dexp( -(2.0d0*ran()-1.0d0) / 2.0d0 ) / dsqrt(2.0d0 * pi)
	  md_v(i,3) = dexp( -(2.0d0*ran()-1.0d0) / 2.0d0 ) / dsqrt(2.0d0 * pi)
	  md_sumv(:) = md_sumv(:) + md_v(i,:)
	end do
	! Remove velocity shift
	md_sumv = md_sumv / md_natoms
	do i=1,md_natoms
	  md_v(i,:) = md_v(i,:) - md_sumv(:)
	end do
	! Rescale velocities for desired temperature
	ke = 0.0d0
	do i=1,md_natoms
	  ke = ke + md_mass(i) * sum(md_v(i,1:3)*md_v(i,1:3))
	end do
	ke = 0.5d0 * ke
	! Calculate instantaneous temperature
	! If ke is in units of [g / mol-1   Angstroms**2 / ps**2] then ke is in units of 10 J mol
 	t_instant = ke / ( 1.5d0 * md_natoms * kb_10jmolk)
	md_v = md_v * sqrt(md_Ttarget / t_instant)
	! Calculate previous atomic positions (for use by basic Verlet algorithm)
	md_rlast = md_r - md_v * md_timestep
	end subroutine md_initialise

	subroutine md_writeframe(newrun)
	implicit none
	integer :: i
	logical :: newrun
	if (newrun) then
	  open(unit=55,file="md.traj.xyz",status="replace")
	else
	  open(unit=55,file="md.traj.xyz",access="append",status="old")
	end if
	write(55,*) md_natoms
	write(55,"(a,es12.5)") "time = ", md_timestep * md_step
	write(55,"(a8,3f12.6)") (md_atomnames(i), md_r(i,1:3), i=1,md_natoms)
	close(55)
	end subroutine md_writeframe

	subroutine md_calc_kinetic()
	implicit none
	integer :: i
	! Determine total kinetic energy and instantaneous temperature of particles
	! As a check, also calculate total velocity of system (md_sumv)
	! KE = SUM(i) [ 0.5 * mass(i) * v(t)**2 ]
	! T = KE / ( 3/2 * N * kb)
	md_ke = 0.0d0
	md_sumv = 0.0d0
	do i=1,md_natoms
	  md_ke = md_ke + md_mass(i) * sum(md_v(i,1:3)*md_v(i,1:3))
	  md_sumv(1:3) = md_sumv(1:3) + md_v(i,1:3)
	end do
	! Finalise total kinetic energy (J mol-1)
	md_ke = md_ke * 0.5d0
	! Determine instantaneous temperature
	! If ke calculated from units of [g / mol-1   Angstroms**2 / ps**2] then ke is in units of 10 J mol
	md_Tcalc = md_ke / ( 1.5d0 * md_natoms * kb_10jmolk)
	! Convert kinetic energy from 10 J/mol to eV
	md_ke = md_ke / 9648.536d0
	end subroutine md_calc_kinetic

	subroutine md_verlet(forces)
	implicit none
	integer :: i
	real*8, intent(in) :: forces(md_natoms,3)
	real*8 :: dt, dtsq, twodt, a(3)
	! Propagate system according to Verlet algorithm
	! Forces are supplied in units of eV/Anstrom**2, so convert to J mol-1 / Angstrom**2 along the way
	dt = md_timestep
	dtsq = dt*dt
	twodt = 2.0d0 * md_timestep
	! r(t+dt) = 2*r(t) - r(t-dt) + f(t)/m * dt**2
	do i=1,md_natoms
	  a(1:3) = (forces(i,1:3) * 9648.536d0) / md_mass(i) 
	  md_temp(i,1) = 2.0d0*md_r(i,1) - md_rlast(i,1) + a(1) * dtsq
	  md_temp(i,2) = 2.0d0*md_r(i,2) - md_rlast(i,2) + a(2) * dtsq
	  md_temp(i,3) = 2.0d0*md_r(i,3) - md_rlast(i,3) + a(3) * dtsq
	end do
	! v(d+dt) = (r(t+dt) - r(t-dt)) / 2dt
	do i=1,md_natoms
	  md_v(i,1) = (md_temp(i,1) - md_rlast(i,1)) / twodt
	  md_v(i,2) = (md_temp(i,2) - md_rlast(i,2)) / twodt
	  md_v(i,3) = (md_temp(i,3) - md_rlast(i,3)) / twodt
	end do
	! Store new arrays and replace old
	md_rlast = md_r
	md_r = md_temp
	! Fold new atomic positions into the unit cell
	call md_foldatoms()
	! PBC old positions so we can calculate velocity properly on next iteration
	call md_pbcatoms(md_rlast, md_r)
	end subroutine md_verlet

	subroutine md_velocity_verlet(forces)
	implicit none
	integer :: i
	real*8, intent(in) :: forces(md_natoms,3)
	real*8 :: dt, dtsq, a(3)
	! Propagate system according to Velocity Verlet algorithm
	dt = md_timestep
	dtsq = dt*dt
	! r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2
	! v(t+dt/2) = v(t) + 0.5*a(t)*dt
	! a(t+dt) = F(t+dt)/m
	! v(t+dt) = v(t) + 0.5*a(t+dt)*dt
	do i=1,md_natoms
	  a(1:3) = forces(i,1:3) / md_mass(i) * 96485.36d0

	  ! ..and finally velocities again (by second half-step)
	  md_v(i,1:3) = md_v(i,1:3) + 0.5d0*a(1:3)*dt

	  ! Propagate positions (by whole step)...
	  md_r(i,1:3) = md_r(i,1:3) + md_v(i,1:3)*dt + 0.5d0*a(1:3)*dtsq
	  ! ...velocities (by half step)...
	  md_v(i,1:3) = md_v(i,1:3) + 0.5d0*a(1:3)*dt

	end do
	! Fold new atomic positions into the unit cell
	call md_foldatoms()
	end subroutine md_velocity_verlet

	subroutine md_foldatoms()
	use cellutils
	implicit none
	integer :: i
	do i=1,md_natoms
	  call fold(md_cell, md_icell, md_r(i,1), md_r(i,2), md_r(i,3))
	end do
	end subroutine md_foldatoms

	subroutine md_pbcatoms(r,reference)
	use cellutils
	implicit none
	real*8, intent(inout) :: r(md_natoms,3)
	real*8, intent(in) :: reference(md_natoms,3)
	real*8 :: delta(3)
	integer :: i
	do i=1,md_natoms
	  call mimvec(md_cell, md_icell, r(i,1), r(i,2), r(i,3), reference(i,1), reference(i,2), reference(i,3), delta(1), delta(2), delta(3))
	  r(i,:) = reference(i,:) + delta(:)
	end do
	end subroutine md_pbcatoms

	end module mdroutines
